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Abstract 

A generalization of the Cucker-Smale model for collective animal behaviour is investigated. The model 
is formulated as a system of delayed stochastic differential equations. It incorporates two additional pro¬ 
cesses which are present in animal decision making, but are often neglected in modelling: (i) stochasticity 
(imperfections) of individual behaviour; and (ii) delayed responses of individuals to signals in their envi¬ 
ronment. Sufficient conditions for flocking for the generalized Cucker-Smale model are derived by using 
a suitable Lyapunov functional. As a byproduct, a new result regarding the asymptotic behaviour of 
delayed geometric Brownian motion is obtained. In the second part of the paper results of systematic 
numerical simulations are presented. They not only illustrate the analytical results, but hint at a some¬ 
how surprising behaviour of the system - namely, that an introduction of intermediate time delay may 
facilitate flocking. 

Keywords: Cucker-Smale system, flocking, asymptotic behaviour, noise, delay, geometric Brownian 
motion. 


1 Introduction 

Collective coordinated motion of autonomous self-propelled agents with self-organization into robust patterns 
appears in many applications ranging from animal herding to the emergence of common languages in primitive 
societies [26] . Apart from its biological and evolutionary relevance, collective phenomena play a prominent role 
in many other scientific disciplines, such as robotics, control theory, economics and social sciences |3 EH H3]. 
In this paper we study the interplay of noise and delay on collective behaviour. We investigate a modification 
of the well known Cucker-Smale model [30 with multiplicative noise and reaction delays. 

We consider N £ N autonomous agents described by their phase-space coordinates (xi(t),Vi(t)) £ R 2d , 
i = 1,2,..., TV, where Xi = Xi{t) £ R d (resp. Vi = Vi(t) £ R d ) are time-dependent position (resp. velocity) 
vectors of the i-th agent. The governing equations are given as the following system of delayed ltd stochastic 
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differential equations 


d Xi = Vidt, (1) 

X N N 

dWi = JyYl ^ & ~ %) dt + y H V'ij fa - «*) dB?, (2) 

j=i j=i 

where the delayed velocity is given by Vipt) = Vi{t — r) and r > 0 is the reaction delay. The parameters A > 0 
and Oi € K, % = 1, 2,..., N, measure the alignment and noise strength, respectively, and dBf, i = 1, 2,..., N, 
are independent d-dimensional white noise vectors. In general, the communication rates ipij are functions of 
the mutual distances |Xi — Xj\, however, in most of our paper we will consider them as given functions of time 
satisfying certain assumptions. The standard Cucker-Smale model mm is a special case of equations © © 
for <jj = 0 and r = 0. Our aim is to investigate equations © © for general values of reaction delay r and 
noise strength parameters er*, i = 1,2 ,,N. 

The Cucker-Smale model was introduced and studied in the seminal papers mm, originally as a model 
for language evolution. Later the interpretation as a model for flocking in animals (birds) prevailed. In 
general, the term flocking refers to the phenomena where autonomous agents reach a consensus based on 
limited environmental information and simple rules. The Cucker-Smale model is a simple relaxation-type 
model that reveals a phase transition depending on the intensity of communication between agents. Using 
Ui — 0 and t = 0 in © ©, we can write the standard Cucker-Smale model as the following system of 
ordinary differential equations 


Xi = Vi, 

A N 

Vi = - Vi), for i= 1,2, 

j =i 


(3) 

(4) 


where the dots denote the time derivatives. We note that the scaling by N~ x in © is significant to obtain 
a Vlasov-type kinetic equation in the mean-field limit N —> oo, see, for example |2Sj . The communication 
rates ipij introduced in El HI and most of the subsequent papers are of the form 


ipij = 1p(\Xi - Xj\) 


with ip{s) 


I 

(i + sY’ 


where (3 > 0. 


(5) 


If j3 < 1/2, then the model exhibits the so-called unconditional flocking, where for every initial configuration 
the velocities Vi(t) converge to the common consensus value N~ 1 y/ , Wj(0) as t —> oo. On the other hand, 
with /3 > 1/2 the flocking is conditional, i.e., the asymptotic behaviour of the system depends on the value 
of A and on the initial configuration. This result was first proved in Eim using tools from graph theory 
(spectral properties of graph Laplacian), and slightly later reproved in [28j by means of elementary calculus. 
Another proof has been provided in m, based on bounding © © by a system of dissipative differential 
inequalities, and, finally, the proof of [J is based on bounding the maximal velocity. 

Various modifications of the generic model © © have been considered. For instance, the case of singular 
communication rates ip(s) = l/s? was studied in HHIH]. Motsch and Tadmor m scaled the communication 
rate between the agents in terms of their relative distance, so that their model does not involve any explicit 
dependence on the number of agents. The dependence of the communication rate on the topological rather 
than metric distance between agents was introduced in m ■ The influence of additive noise in individual 
velocity measurements was studied in m and |29] . Stochastic flocking dynamics with multiplicative white 
noises was considered in |lj. Delays in information processing were considered in |l8j . however, their analysis 
only applies to the Motsch-Tadmor variant of the model. 

In this paper, we are interested in studying the combined influence of noise and delays on the asymptotic 
behaviour of the Cucker-Smale system. In particular, we derive a sufficient condition in terms of noise 
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intensities <j.; and delay length r that guarantees flocking. Our analysis is based on a construction of a 
Lyapunov functional and an estimate of its decay rate. To prove our main results, we make an additional 
structural assumption about the matrix of communication rates which, loosely speaking, means that the 
communication between agents is strong enough. 

The paper is organized as follows: In Section[2]we show how our model JT]) © is derived from the Cucker- 
Smale model ([3]) (]4|) and define what is meant by flocking. Moreover, we consider a simplified version of the 
model to provide an intuitive understanding of what qualitative properties may be expected. In Section [3] 
we derive a sufficient condition for flocking in terms of the parameters A, cq and r, based on a micro-macro 
decomposition and construction of a Lyapunov functional. Moreover, as a byproduct of our analysis, we 
provide a new result about the asymptotic behaviour of delayed geometric Brownian motion. Section [I] is 
devoted to a systematic numerical study of the model. First, we focus on simulation of delayed geometric 
Brownian motion, in particular, we study the dependence of its asymptotic behaviour on the delay and noise 
levels. Then, we perform the same study for system m ©. This leads to the interesting observation that, 
for weak coupling and small noise levels, an introduction of intermediate delays may facilitate flocking. A 
systematic study of this effect concludes the paper. 

2 The stochastic Cucker-Smale model with delay 

In order to make the generic model d3| d3| more realistic, we amend it with two additional features. First, we 
note that measurements in the real world are subject to errors and imprecisions that are typically modeled 
in terms of white noise. In particular, we assume that the state (velocity) of agent j measured by agent i is 
given by the expression 


Wi;j = Vj + Ki(vj - Vi) dB\, for i, j = 1, 2,..., TV, (6) 

where re* > 0 represents the imprecision of i 's measurement device, and B\ are independent identically 
distributed d-dimensional Brownian motions with zero mean and the covariance relations 

E [Bj ■ B*j] = d 5ijSt s t, for i,j = 1,2,. .., N and t, s > 0, 

with Sij = 1 iff i = j and Sij = 0 otherwise, and similarly for <5 ts . 

Note that the multiplicative structure of the noise term ensures that uii-i = v t . Substituting uji-j given by 
© for Vj in 0 and defining eq := A Hi, we obtain the following system of stochastic differential equations 
(SDEs) for velocities 


A n n 

Avi = + ~ v i) dB h for i = l,2,...,N, 

i =i t=i 

with A a positive coupling strength. 

The second amendment is the introduction of delays, motivated by the fact that agents react to information 
received from their surroundings with some time lag. However, we assume that information propagates 
instantaneously, so the delay does not depend on the physical distance between agents. For simplicity, we 
assume the reaction lag to be the same for all agents, so that at time t they react to information perceived 
at time t — t for a fixed r > 0. 

Convention 1 Throughout the paper, we denote by Vi the quantity Vi evaluated at timet, i.e., Vi = Vj(t), and 
by Vi the same quantity evaluated at time t — T, i.e., Vi = Vi(t — r). We will also write x = (x\, X 2 ,..., Xn) € 
( reS p. v = (vi,V 2 , ■ ■ ■ ,vn ) G R dxAr ) for the vectors of locations (resp. velocities ) of the agents. 
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In general, the communication rates ipij > 0 may be functions of the mutual distances |Xi — Xj |. However, 
our analysis is based on a certain structural assumption about the communication matrix an d the 

particular form of the dependence on the mutual distances is irrelevant. Therefore, we consider the rates 
ipij = as given adapted stochastic processes, so that © decouples from ©. Moreover, we assume that 

ipij are uniformly bounded, 

0 < < 1, for t > 0, i, j = 1, 2,..., N, almost surely. (7) 

Thus, we finally arrive at the stochastic system of delayed differential equations that we will study, 


N 


N 


dv i = V’y { v j - Vi) dB \> 

J=1 3 =1 

which is supplemented with the deterministic constant initial datum v° £ R dxJV , 

Vi(t) = v°, for f e (-r, 0], i = 1,2,..., N. 

Let us note that we interpret the noise term in © in terms of the Ito calculus I2i nn. 


( 8 ) 


(9) 


Theorem 1 The stochastic delay differential system © with initial condition © admits a unique global 


solution v = v(£) on [—r, oo) which is an adapted process with E J_ r |v(t)| 2 d< 
martingale. 


< oo for all T < oo, i.e., a 


Proof: The proof follows directly from Theorem 3.1 and the subsequent remark on p. 157 of m ■ Indeed, 
© is of the form 

dv(f) = F{t, v(t — r)) d t + G(t, v(t — r)) d B l 

for suitable functions F and G. In particular, the right-hand side is independent of the present state v(£), 
so that the solution can be constructed by the method of steps. The second order moment is bounded on 
(—r, T) because of the linear growth of the right-hand side of © in v. 


We now define the property of asymptotic flocking for the solutions of © ©• 

Definition 1 We say that the system © exhibits asymptotic flocking if the solution fv{t))t>o for any initial 
condition © satisfies 


lim 

t—foo 


E[ui(£)] -E[i>j(f)] 


= 0 


for all i, j = 1, 2,..., N, 


where E[-] denotes the expected value of a stochastic process. 

2.1 Simplified case with ^ = 1 

To get an intuition of what qualitative properties we may expect from the solutions of ©, we consider the 
case when the communication rate is constant, i.e., ipij = 1; in other words, we set f = 0 in ©• We also 
assume that eq is equal to the same constant a £ R for all z = 1,2,..., TV, i.e. cq = a, and, moreover, that 
vi = v° for some v° £ and all i = 1, 2,..., TV. Then, defining V c (t) := jj Vi(t), we obtain 


N 


i= 1 
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Since, by assumption, V c (t) — Vi(t) = 0 for t € (—r, 0], we have V c (t) = v° for all t > 0. Consequently, © 
decouples into N copies of the delayed SDE 

dw =—\w dt — cr w dB f , (10) 

where we denoted w := Vi — v° for any i = 1,2,..., TV. We are not aware of any results concerning the 
asymptotic behaviour of m- The method developed in [3] suggests that 

r°° i 

lim E riw;(t)| 2 l =0 if and only if / r\(t) 2 dt < —, 

t-i-oo J Q cH 

where r\ is the fundamental solution of the delayed ODE 

w = —A w, (11) 

i.e., formally, r\ solves ED subject to the initial condition w(t ) = f° r t £ (—t, 0]. The fundamental 

solution r\ can be constructed by the method of steps j25], however, evaluation of its L 2 (0, oo)-norm is an 
open problem. From this point of view, the analysis carried out in Section [3] provides new and valuable 
information about the asymptotics of (flT)l) , see Section 13.41 Let us note that setting r = 0 in the above 
criterion recovers the well-known result about geometric Brownian motion : the mean square fluctuation 
E[|ui(f)| 2 ] tends to zero if and only if a 2 < 2A. 

Finally, for the convenience of the reader, we give an overview of the qualitative behaviour of solutions 
to ED with A > 0, subject to a constant nonzero initial datum (see, e.g., Chapter 2 of [25]): 

• If At < 1/e, the solution monotonically converges to zero as t —> oo, hence no oscillations occur. 

• If 1/e < At < 7t/ 2, oscillations appear, however, with asymptotically vanishing amplitude. 

• If At = 7t/ 2, periodic solutions exist. 

• If At > 7t/ 2, the amplitude of the oscillations diverges as t —> oo. 

Hence, we conclude that the (over)simplified model (fill) , corresponding to the delayed Cucker-Smale system 
with ip = 1 and no noise, exhibits flocking if and only if At < n/2. In the next Section we derive a sufficient 
condition for flocking for the general model ©. 

3 Sufficient condition for flocking 

In this section we derive a sufficient condition for flocking in © according to Definition [Q Our analysis 
will be based on a construction of a Lyapunov functional that will imply decay of velocity fluctuations for 
suitable parameter values. However, we will have to adopt an additional structural assumption on the matrix 
of communication rates {ipij)^- 1 - 

Before we proceed, let us shortly point out the mathematical difficulties that arise due to the introduction 
of delay and noise into the Cucker-Smale system. The “traditional” proofs of flocking of model © ©, for 
instance mmmm, rely on the monotone decay of the kinetic energy (velocity fluctations) of the form 

N . N N 

-j t J2 n 2 = ^\ v i - «j-i 2 < o. 

i— 1 i= 1 j= 1 

However, this approach fails if processing delays are introduced, since for © without noise (i.e., all (7^ = 0), 
we have 

, n . n N 

77 n 2 = • (% - ^)- 

i=l i= 1 j=1 
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One then expects the product (vi — vf) ■ (u, — vf) to be nonnegative for r > 0 small enough, however, it is 
not clear how to prove this hypothesis. 

The introduction of noise leads to additional difficulties - in particular, the classical bootstrapping argu¬ 
ment 00E2 for fluctuations in velocity fails in this case. Similarly as in m , we circumvent this problem 
by adopting, in addition to the boundedness m, a structural assumption about the matrix of communication 
rates. We define the Laplacian matrix A(f) £ S> NxN by 

Ayj .— ipiji for i j, Ayy -— ^ ? (12) 

and note that A is symmetric, diagonally dominant with non-negative diagonal entries, thus it is positive 
semidefinite and has real nonnegative eigenvalues. Due to its Laplacian structure, its smallest eigenvalue is 
zero [6, . Let us denote its second smallest eigenvalue (the Fiedler number) p .2 (t). Our structural assumption 
is that there exists an £ > 0 such that 

H 2 (t) > £ > 0, for t > 0, almost surely. (13) 

This can be guaranteed for instance by assuming that the communication rates are uniformly bounded away 
from zero, ipij(t) > ip > 0 , since there exists a constant c > 0 such that /^(i) > cip, see Proposition 2 in (6j. 

Moreover, we assume that the matrix of communication rates is uniformly Lipschitz continuous in the 
Frobenius norm, in particular, there exists a constant L > 0 such that 

||A(t) — A(i — t)|| f < Lt for t > 0, (14) 

where ||-|| f denotes the Frobenius matrix norm. 

To ease the notation and without loss of generality, we will consider the one-dimensional setting d = 1, 
i.e., Vi(t ) £ R and v(f) £ R w , where N is the number of agents. Then, with the definition (fT2l) . we put (0 
into the form 

d^i =--^r(Av)jdt + ^(Av)jdB|, i = 1,2,..., AT. (15) 

Our main result is the following. 


Theorem 2 Let A be given by m satisfying o, m and (ED- Let the parameters A > 0 and cr^x := 
max{cr^, (72,..., cf[} satisfy 

^max < A, (16) 


then there exists a critical delay t c = r c (A, cr max , L, £) > 0, independent of N, such that for every 0 < r < r c 
the system (1151) exhibits flocking in the sense of Definition [H 

Moreover, if the matrix of communication rates A is constant, i.e. ED holds with L = 0, then t c is of 
the form 


t c 


1 

W 


— (7 


2 

max 



+ I2 (A 


- ( 7 2 ) 2 

max/ 


(17) 


Remark 1 The system (1151) with constant communication matrix A can be seen as a linearization of the 
system ©-© about the equilibrium Vi = Vo for i = 1, 2,..., TV with some vq £ R. Note that in this case the 
formula El for the critical delay t c does not depend on the particular value of £ in m- 
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3.1 Micro-macro decomposition 

We introduce a micro-macro decomposition [28], na which splits m into two parts: macroscopic, that 
describes the coarse-scale dynamics, and microscopic, that describes the fine-scale dynamics. The macroscopic 
part for the solution is set to be the mean velocity V c (t), 

1 N 

Vc(t) :=( 18 ) 

V 2=1 

The microscopic variables are then taken as the fluctuations around their mean values, 

Wi(t) := Vi(t) - V c (t), for i = 1, 2,..., N. (19) 

We denote w(t) = (wi,W 2 , ■ ■ ■, wn ) £ M. N . Then we have 

w(t) = v(t) — V c (t) e, where e := (1,1,..., 1) T £ R w . (20) 

Since e is the eigenvector of A corresponding to the zero eigenvalue, we have Aw = Av. Then (1151) can be 
rewritten as follows 

dwi = ~(Aw), dt + ^(Aw)j d B\ - dV c . (21) 

The macroscopic variable V c satisfies the following lemma. 


Lemma 1 Let (y(t))t>o be a solution of (J8[) subject to the deterministic constant initial datum ©. Then 


E[V),(t)] = V),(0) for t > 0 and E J_ t |V' c (t)| 2 dt 


< oo for all T < oo. 


Proof: The boundedness of E 
property of v(t) provided by Theorem [l] Using (|T2|) . we have 


f_ r |U c (i)| 2 dt follows directly from the definition (fTSl) and the martingale 


N 

^(Av)j = e T Av = 0. 

i =1 

Summing equations (fl5l) . i = 1,2,..., N, using (TTH1) and Aw = Av, we obtain that the macroscopic dynamics 
is governed by the system 

1 N 

dVc = H ai ( AW ^ dB i ' ( 22 ) 

i -1 

After integration in time this implies 


1 N 

nv c (t)] = v, :( 0 ) + F ^ 

i=l 


; E 


[ (Ms - t)w(s - r))i d B- 
Jq 


for t > 0. 


Since f(s) := (A(s — r)w(s — r))j is a martingale, we have E 
[T5]). Thus we obtain E[U c (t)] = V c (0). 


fof(s)dBf 


0 (see Theorem 5.8 on p. 22 of 


Remark 2 Note that m and (|22l) are expressed in terms of the xv-variables only and so they form a closed 
system, which is equivalent to m- 
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Clearly, due to m, we have w • e = EE w i — 0- Consequently, it is natural to introduce the decompo¬ 
sition R w = (e) © (e)- 1 -, where e is given by (EU1) . We then have w (t) £ (e)* L for all t > 0. 

Lemma 2 Let A £ R NxN , N >2, be the matrix defined in m and assume that (ED and JED hold. Then: 

(a) The maximal eigenvalue of A is bounded by 2 ( N — 1). 

(b) We have |Au| 2 < 2 ( N — 1) u T Au for any vector u £ R w . 

(c) We have £ |w| 2 < w T Aw < 2 (N — 1) |w| 2 for any vector w £ (e) 2 -. 

(d) For any vectors u, w £ (e) 1 - and S > 0 we have 

u t Aw < — u t Au H— w t Aw. 

~ 26 2 


Proof: 


(a) The claim follows from the Gershgorin circle theorem. Indeed, since 0 < ipij < 1, the diagonal entries 
satisfy 0 < An < N — 1, and l A y| = A i* for all i = 1,2,..., AT. 


(b) The smallest eigenvalue of A is zero with the corresponding eigenvector e. The second smallest eigenvalue 
fj ,2 (Fiedler number) is assumed to be positive by (1131) . Thus, A is a symmetric, positive operator on the 
space (e)- L and there exists an orthonormal basis of (e) 2 - composed of eigenvectors £ 2 ,£ 3 , ■ ■ • ,Zn of A 
corresponding to the positive eigenvalues M 2 ,M3, ■ • •, mjv- Then, every vector u £ R^ can be decomposed 
as 


u = 


(u • e) e 


N 


+ 7 X u ■ ii) Zi- 


i -2 


(23) 


Thus, due to the above bound on the eigenvalues 0 < Mi < 2 (N — 1), we have 


N 


N 


= ‘ &)W - 2 ( N - X ) ‘ = 2 ( N - 1) U 3 


Au. 


i =2 


i =2 


(c) If w G (e)- 1 , then (f23l) implies 

N N 

’ w = ^2( w 'Zi)Zi, and W T Aw = E( w '£i) 2 ^- 

i =2 i =2 

Since nonzero eigenvalues are bounded from below by £ (using (fEll l and from above by part (a) of this 
lemma, we obtain 

N N N 

£ |w| 2 = £ 5> • Zi) 2 < • ZiftH < 2 {N - 1) £(w • Zi) 2 = 2(N-1) |w| 2 . 

i —2 i —2 i—2 

(d) With the orthonormality of the basis Z 21 Z 3 , ■ ■ ■ ■ Zn an d the positivity of the eigenvalues M 2 , M 3 , • • •, Hn, 
we have by the Cauchy-Schwartz inequality 

( N \ T /N \ N 

E^ U ' &) ) A ( ) = 

i =2 / \i =2 / i =2 

( N \i/N \ 5 

5Z(u-^) 2 /i*J (E( W -^J = (u T Au)5(w T Aw)5, 
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and with any <5 > 0, 


(u 1 Au) 2 (w t Aw) ^ < -L u t Au + ^ w t Aw. 

2 o 2 


3.2 Lyapunov functional 

The proof of Theorem[2]relies on estimating the decay rate of the following Lyapunov functional for (1211) (l22l) . 

Jf(t) := |w(t)| 2 + -^2 £ |A(s)w(s)| 2 ds (24) 

+ ^2 J t J g |A(s-r)w(s-r)| 2 ds, 

where p, q are positive constants depending on A, r and eq. 


Lemma 3 Let the assumptions of Theorem [2] be satisfied. Then there exist positive constants p, q and e such 
that for every solution (w(f)) t >o of (1211) - (|22l) the Lyapunov functional (l24l) satisfies 


_d 
d t 


< --e 


w t Aw 


(25) 


Proof: We apply the Ito formula to calculate d Wi(t) 2 . Note that the Ito formula holds in its usual form also 
for systems of delayed SDE, see page 32 in [12] and 1130 EDI- Therefore, we obtain 


d Wi(t) 2 = Wi(Aw)i d t + jti)i(Aw)i d B\ - 2wi dV c 

2 ^ N N 

+ ^2^)i dt+ fy4 EE Vj°k(&w)j(Aw)k d B) d B[ 

3 =1 k=l 

-2 (^(AfiOidBj) dV c . 

With the identity dd Bf. = Sjk dt (formula (6.11) on p. 36 of [19]), we have 


1 

TV 1 


N N 

E y (TjCTk{Aw)j(Aw)k d Bj d B* k 

3=1 fe = 1 


1 


N 

3=1 


) 2 df. 


Consequently, summing over i. using w • e and the identity 


E d B i) dVc= ^y dt, 


d|w(f)| 2 = f — ^w r Aw + 


TV — 1 
N 3 


N 

E 

2=1 


N 


G i (Aw) i d£ + — > (TiWi(Aw)i dBj . 


N ^ 

2 = 1 


we obtain 
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Consequently, we have 


2A ~ N — 1 ~ \ 2 ~ 

dJf(t) = ( - —w t Aw + —^ 3 - ^ erf (Aw)? ) di + — 53 ViWi(A-w)i d B\ 


N 3 ' N- 

2=1 / 2=1 


+ 


TV 2 


(|Aw| 2 - |Aw| 2 ) dt + ^ ^r|Aw| 2 - J* |Aw(s - r)| 2 dsj df. 


(26) 


Our goal is to estimate -gyEfJz?(f)] from above. First of all, we note that by the elementary property of the 
Ito integral (Theorem 5.8 on p. 22 of [T9]), 


E 


1 N 

— 53°‘i w i( Aw )idS- 


= 0. 


For the first term of the right-hand side in (1261) . we write 


2 A 


2 A 


2 A 


— -—w J Aw =- w 1 Aw H-w 2 A(w — w) 

TV TV TV y ’ 


and apply Lemma H)[d) with 5 > 0, 
2 A 


r. ~ A^" 1 

—w" A(w - w) < w 


t Aw + -^-(w — w) t A(w — w). 


Using Lemma|2]c), we have 


A <5 


TV- 1 


(w — w) t A(w — w) < 2 A <5 

TV 1 ' v ’ “ TV 

Now we write for w — w, componentwise, using (1211) . 


w — w . 


(27) 


= 


/ dtUi(s) 

J t—T 


= -4 / (A(s — t)w(s — r)), ds + ^ 


TV 


TV 


(s - t)w(s - r))jdS? 


dU c (s). 


Thus, we have for the expectation of the square 

A 


E|wj - ■ 


->l\ JS 


< 3E 


TV 


+3E 


s — t)w(s — r))j ds 

dU c (s) 


't — T 

■t -.2 


+ 3E 


-12 


^ (A(s-T)w(s-T))jdB? 


_*/ t — T 


An application of the Cauchy-Schwartz inequality and Fubini’s theorem for the first term of the right-hand 
side yields 


3E 


4 h ( A ( s — r ) w (s — r))i ds <44 r j E[|(A(s-r)w(s-r)) i | 2 ] ds. 
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For the second term we use the fundamental property of the Ito integral (Theorem 5.8 on p. 22 of [T9ll. 


3E 


^It ( A ( s _ T ) w ( s _ ^ djB ’ ? = ~^tj t E [K A ( s ^ T ) w ( s ^ r )) i l 2 ] ds - 


Similarly, the third term is estimated as 

-,2 


3E 


dV' c (s) 


N 4 


E 


J2 / a j( A (s - r)w{s ~ t))j dB^ 
j= i Jt ~ T 


t—T 

t 


< 


< 


S I 

J t Vj(Hs - t)w(s - T))j dB. 

J2 [ E[|(A(s-r)w(s-r)) i | 2 ] ds. 

A — 1 J t — T 


.2 

max 


Thus, we get from (12711 . estimating A_J- < 1, 


XS ^ 

— E 

N 


n — 
6A S 


N rt 


(w - w) t A(w - w) < - 2 - (A 2 r + 2cr^ ax ) / E[|(A(s - r)w(s - r))i| 2 ] ds. 

*=l • 7t - T 

An application of Lemma [2Jb) gives 


Q TTMI A_l2l ^ 2 9 (A” 1 ) njir T 


N2 n^]< N2 


E[w J Aw] < ^E[w t Aw], 


To balance this term with — 2Aw T Aw, we use assumption (1141) and Lemma [2](c) in 


f T Aw = w t (A — A)w + w t Aw 


< Lrw T w + w r Aw < (if 1 r + l) w t Aw. 

Collecting all the terms in (OBI) finally leads to 

-^E[_S?(t)] < i_ [_2A +Ad" 1 +2g(ir 1 r + l)] E[w T Aw] + ^ (cr^ ax + pr- 9 )E|A 


N 


N 2 


We set 


+ ^2 ( 6 A<5(A 2 r + 2c r ^ ax )-p)^' E|A(s-r)w(s-r)| 2 ds. 


p = 6 XS (A 2 r + 2cr^ ax ), g = cr^ ax +pr, 


then the above expression simplifies to 


E[JSf(i)] < -1 [-2 A + A (T 1 + 2 9 (irV + l)] E[w T 


(28) 


(29) 


We want —2 A + Ad -1 +2 q ( L£~ 1 t + l) <0. Substituting (l28l) into this inequality leads to a third order 
polynomial inequality in r. This polynomial has all positive coefficients but the zero order one, which is 
c 0 := 2cr 2 lax + d _1 A — 2A. If (fl6l) is satisfied, then choosing S > 0 such that 


S < 2A (A - <r max ) 
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makes Cq negative. Consequently, there exists a r c > 0 such that for any 0 < r < r c , 

-2X + XS- 1 +2q(L£~ 1 T+l) =-e<0. 


This completes the proof of (1251) . 

It remains to study the case when A is a constant matrix, i.e. L = 0 in (fldl) . Then (1291) simplifies to 
■^■E[2z?(£)] < -jy (—2 A + A<5 _1 + 2 q) E[w t Aw] 

and we have to find r such that —2 A + A S~ x + 2 q < 0. Again, substituting (1281) for p and q leads to 

_ X(2-S- 1 )-2al ax 

12 A 5 (A 2 t + 2cr^ ax )' 


(30) 


The maximum value of the right hand side is obtained for 5 = A(A — o’max) 1 which is positive because of 
the first inequality in (flTl) . Substituting S = A(A — cr^^) -1 into (l30l) . we obtain 


Finally, resolving in r leads to 


T < 


(A-<4 ax ) 2 


12 A 2 (A 2 r + cr^ ax )' 


r < 


1 


— a 


2 

max 






If the above sharp inequality is satisfied, there exists an £ > 0 such that —£ = —2A + A 5 1 + 2 q and, 
consequently, (B5l) holds. 


3.3 Proof of Theorem [2] 

An integration of (1251) in time gives 

r* r i 

E[|w| 2 (f)] < E[2Sf(t)] = E[2Sf(0)] + / — E[2S?(s)]ds 

Jo dt 

< E[«£?(0)] -Jjj E[w t (s)A(s)w(s)] ds. 

An application of Lemma [He) gives then 

E[|w| 2 (t)] < E[2Sf(0)] - ^ f E[|w(s)| 2 ] ds, 

so that the last integral is convergent as t —> oo and, consequently, lim^oo E[iUj(i)] = 0 for all i = 1, 2,..., N. 
Using (HU), we obtain 


lim |E[uj(£)] — E[uj(i)]| = lim |E[wj(f)] — E[i«j(t)]| = 0, 


t—t OO 


t—> OO 


and we conclude that asymptotic flocking in the sense of Definition [T| takes place. □ 
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Remark 3 The original definition of asymptotic flocking 00 involves, in addition to our Definition^ also 
the group formation property for © © given by 

sup |E[xj(t)] — E[X c (t)]| < oo, for all i = l,2,...,N, 

t> o 


where X c {t) 


1 sr^N 

N 1 


Xi(t). 


The standard way 00121 na of proving this result would be to estimate 


|E[*i(t)]-E[X c (t)]|< |®?-x c (0) 


E[|u>j(s)|] ds 


and employ a bootstrapping argument to show that / Q °° E[|wi(s)|] ds < oo. However, as noted above, it is not 
clear how to apply the bootstrapping argument in our setting. Note that we have J^° E[|u?i(s)| 2 ] ds < oo, but 
it does not imply / 0 °° E[|iyj(s)|] ds < oo. 


Remark 4 Let us note that the flocking conditions (I16H and 03. expressed in terms of A and <r max , translate 
into 


^ K Lx < 1| T < ^max + Y K iax + J^2 _ ^ K max) 2 

in terms of A and re max := cr max /A, where K max can also be written as K max = max{K 1; k 2) ■ • •, kjv} for n it 
i = 1, 2,..., N, defined by ©. 


3.4 Application to asymptotic behaviour of delayed geometric Brownian motion 

Our analysis provides information about the asymptotic behaviour of the delayed geometric Brownian mo¬ 
tion Itinil which is, to our best knowledge, new. We just modify the proof of Lemma 0 with the obvious 
simplifications due to the fact that Aft) = 1. This leads to a slight improvement in the flocking condition. 

Lemma 4 Let the parameters A > 0, a G R and r > 0 satisfy 

a 2 < 2A, t < (-2tr 2 + V4^T2(2A^2)2) . (31) 

Then the solutions of the delayed geometric Brownian motion equation m satisfy 

lim E[H 2 ] = 0. 

t—y OO 

Let us note that the above result is suboptimal for the deterministic case. Indeed, setting a := o, m 
reduces to At < \/2/2. However, it is known [25] that solutions of the delayed ODE w = —A w asymptotically 
converge to zero if At < n/2. On the other hand, if there is no delay, i.e., t = 0, the condition (lUTT) reduces 
to a 2 < 2A, which is the sharp condition for asymptotic vanishing of mean square fluctuations of geometric 
Brownian motion. 
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4 Numerical experiments 

We provide results of numerical experiments for the models considered in this paper with focus on their 
asymptotic behaviour. First, we illustrate that one has to be cautious when interpreting the numerical 
results as indications about the “true” asymptotic behaviour of the solution, because implementations of 
Monte-Carlo algorithms for geometric Brownian motion lead to systematic underestimation of the moments 
of the true solution, see Section 14.11 Keeping this systematic defect in mind, we will resort to weak methods 
for simulation of our SDEs and study their numerical asymptotic behaviour. In Section [4. 2 1 we resort to the 
delayed geometric Brownian motion (flOl) . which can be seen as a toy model of CD ©< and find combinations 
of parameter values that guarantee numerical asymptotic decay of the solution. In Section 14.31 we then 
perform numerical simulations of the velocity alignment system (0 with fixed communication rates, and, 
finally, in Section 14.41 we focus on the full system CD © ■ 


4.1 Analysis of the Monte-Carlo method for geometric Brownian motion 

In this section we estimate the systematic error produced by numerical implementations of the Monte-Carlo 
algorithm for geometric Brownian motion without delay. We show that computer simulations underestimate 
the mean square fluctuations of the process due to the fact that the numerical implementation does not 
capture large deviations (extreme outliers), and the error grows exponentially in time. 

Let us consider the one-dimensional Brownian motion with drift, 

dz = (A — ct 2 /2) dt + odB*, z(0) = 0. (32) 

Then, defining v{t) := exp[z(£)], we have by the Ito formula 

dv = Xvdt-\-crvdB t , u(0) = 1. (33) 

For simplicity, we perform a model calculation with 2 A = a 2 > 0, so that 

dz = V / 2AdI? t , du = Xvdt + V^XvdB*. (34) 

Then, the density u(t,x ) of the process z is given by 

Moreover, we have 

E[z(f)] = 0, E[z 2 (f)] = 2 At, E[u(t)] = exp(Af), E[v 2 (f)] = exp(4At). 


We assume that our numerical scheme produces approximations z of the process z that excludes the extreme 
outliers, i.e., Prob(|z(f)| > a(t)) = 0 for some a = a(t). In particular, we consider a properly scaled cut-off 
of the density u(t,x) such that the probability of the extreme outliers Prob(|z(f)| > a(t)) remains constant 
in time. This leads to a(t) = rj a/4 A t for some r) > 0, since 


pOO 

Prob(|z(t)| > r/VdXt) = 2 / u(t,x)dx = erfc( 77 ) 

J rj a/4 A t 


where erfc( 77 ) = exp(— x 2 ) dx is the complementary error function. Consequently, we turn u(t, x) into 

the truncated probability density 


X[— a(t),a(t)] ( x ) 
erf ( 77 ) 


u(t,x ), 


u v (t,x) : 







On Cucker-Smale model with noise and delay 


15 


where X[-a(t),a(t)l is the characteristic function of interval [— a(£),a(£)] and erf(??) = 1 — erfc(^) is the error 
function. Let us denote by z(t) the process with the density u v (t, x) for a fixed 77 > 0, and v(t) := exp (z(t)). 
A simple calculation then reveals that 

/ OO 

exp( 2 a?) du v {x) = 

-OO 

Consequently, since E[u 2 (£)] = exp(4A£), the numerical method produces the relative error 


SgM («*(—„ -y/m)- erfc(, - v^t)) . 


E[u 2 (£)] erfc(— 77 — \/4A £) — erfc(r/ — V4A £) 

E[u 2 (£)] ~ 2erf(ij) ' 

This ratio is equal to one for £ = 0. Using the mean value theorem, we obtain the asymptotic behaviour of 
the ratio for large times, 

E[h 2 (£)l 2 77 

IpWj “ TTrfM “Pi- 4 "). t ^ oc. 

Consequently, any implementation of the Monte Carlo method excluding large deviations will underestimate 
the true value of E[n 2 (£)] by an exponentially growing factor in time. Let us note that this is also true for 
any moment of v and with general parameters A and <7. We illustrate this fact using a numerical simulation. 
We perform a Monte Carlo simulation in Matlab with 10 6 paths of the process (15121) on the time interval 
[0,T] with T = 30 and 10 3 equidistant sampling points. We impose the initial condition z(0) = 0 and the 
parameter values A = 0.5 and <7 = 1. Consequently, z(t) is the Wiener process £?*, and for its numerical 
approximation 5 we use the built-in Matlab procedure normrnd that generates normally distributed random 
numbers. We calculate v(t) := exp (z(t)) and evaluate the mean squared fluctuations E[u 2 (£)]. We plot its 
logarithm as the solid curve in Figure [Ha), compared to the analytical curve log(E[u 2 (£)]) = 4A£ (dashed 
line). We observe the exponential in time divergence of the two curves. This is well described by our formula 
(ED, as illustrated in Figure [lib). For the calculation of the cut-off parameter 77 we use the maximal value 
attained by the actual numerical realization of the stochastic process, i.e., we set 77 := max t 6 ( 0 ,r] ■ We 
then plot the logarithm of the ratio E[u 2 (£)]/E[t; 2 (£)] and the theoretically calculated curve given by the 
right-hand side of (l35l) . We observe a good match between the two curves. 

This systematic discrepancy between the analytical formulas and Monte Carlo simulations originates in 
the heavy tailed distribution of the geometric Brownian motion, and is a well studied topic, see, e.g., the 
survey m ■ Importance sampling and rare-event simulation techniques would be the methods of choice to 
overcome this problem, however, their implementation is beyond the scope of our paper. Instead, we argue 
that the motivation for studying the stochastic system m © and its simplification m comes from the fact 
that they ought to represent models of some real (physical or biological) phenomena. In real life situations 
these extreme events with exponentially low probabilities can be unphysical and the presented Monte Carlo 
simulation can be considered a more appropriate description of reality than the SDEs. We adopt this point 
of view for the forthcoming numerical studies in Sections 14.2114.31 and 14.41 and accept the fact their results 
may not be directly related to statements of Theorem [2] and Lemma [4] 


4.2 Numerical study of delayed geometric Brownian motion 

Using A = 1, the delayed SDE (fllil) can be equivalently written as 

die = — w d£ + cr w d B l , (36) 

where o > 0 and r > 0 are nonnegative parameters. We perform a systematic numerical study of the 
delayed SDE (l36ll to characterize the asymptotic behaviour of its solutions in dependence on the values of 
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Figure 1: (a) Logarithm of the simulated mean squared fluctuations log(E[ti 2 (t)]) (solid line), and the analytical 
result log(E[u 2 (t)]) = 4 A t (dashed line), (b) Logarithm of the ratio of simulated and analytically calculated 
mean squared fluctuations E[u 2 (t)]/E[u 2 (t)] (solid line), and the theoretically calculated curve (dashed line) 
given by the right-hand side of (GU). The Monte-Carlo simulation for z(t) was performed with 10 6 paths of 
the process El) with z(O) = 0, A = 0.5 and a = 1 on the time interval [0, 30] divided into 10 3 equidistant 
sampling points. 


the parameters a and r. In particular, we divide the domain [0, 2] x [0, 2] for (<r, r) into 100 x 100 equidistant 
(cr, r)-pairs. For each pair of the parameter values we perform a Monte Carlo simulation for (1361) with 
Q = 100 paths over the time interval [0,T] with T = 30 and timestep At = 10 -3 . We impose the constant 
deterministic initial condition w(t ) = 1 for t £ (—r, 0]. For discretization of (1351) we use the Euler-Maruyama 
method, i.e., the discrete scheme is 

wt k+ i = w tk -/ktw tk - T + a VAtAfoy, k = 1,2,..., K, (37) 

subject to the initial condition wt = 1 for t < 0. Here K = T / At denotes the total number of timesteps, 
ifc = kAt, and J\fo,i a normally distributed random variable with zero mean and unit variance. Note that the 
values of r are chosen to be integer multiples of At, so that tk — t = ti for some l £ Z. For each (<r, r)-pair 
and each path q of the Monte Carlo simulation we calculate the “indicator” 

q / T/At 

J ^ := n£ Kl 2 

V q—l \ k=(T—l)/At 

where w( is the g-th path in the Monte Carlo simulation of (1371) . The background colour in Figure [2] encodes 
the logarithm of I aT . To define a region of “numerical convergence”, we choose a threshold 0 such that 
I 0 Tc = 0 for the delay r c = ^ that is critical for the problem without noise (cr = 0). In our case this led to 
0 ~ 10~ 2 . The region of “numerical convergence” is marked dark blue in Figure [2] We observe the decrease 
of the critical value of the delay with increasing level of noise. For comparison, the critical values of r given 
by formula (1311) as a function of a are indicated by the solid line. 
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Figure 2: Results of Monte Carlo simulations of the delayed SDE (1561) with Q = 100 paths for ( a , r) £ 
[0,2] x [0,2], The background colour encodes log(/ CT!r ). The region of “numerical convergence” is dark blue. 
The solid line indicates the critical values of r given by formula (l3lT) with A = 1 as a function of a. 


4.3 Numerical study of the system ([ 8 ]) with fixed communication matrix 

We present results of numerical simulations of the system ([5]) in the one-dimensional setting d = 1, where 
we fix the communication rates to ipij = 1 for all i,j = 1,2,... ,7V, i.e., every agent communicates with all 
others at the same rate. Consequently, the communication matrix A has the off-diagonal entries A t j = — 1, 
i ^ j, and A„; = TV — 1. It only has two eigenvalues, 0 and TV. Consequently, its Fiedler number is /X 2 = TV 
and we can choose l := TV in (fT51) . In this setting, we can directly compare our analytical result, Theorem [2l 
with numerical simulations. 

We will be considering even number of agents TV = 2 AT, in particular, TV £ {2,20}, and prescribe the 
initial datum 


Viit) = 



for i = 1, 2,..., K; 

for i = K + 1, K + 2,..., TV; 


(38) 


for t £ (—r, 0]. Although the asymptotic behaviour of the solutions in general depends on the particular choice 
of the initial datum, a systematic study of this dependence is beyond the scope of this paper. Therefore we 
only consider the “generic” choice of initial conditions (l38l) . 

We perform Monte Carlo simulations of the system ([5]) with TV £ {2, 20}, cq = a for all i = 1,2,..., TV 
and A = 1 (other values of A can be achieved by rescaling of a and time). We divide the domain [0,2] x [0,2] 
for (cr, r) into 50 x 50 equidistant (cr, T)-pairs. For each pair of the parameter values we perform a Monte 
Carlo simulation with Q = 100 paths over the time interval [0,T] with T = 30. We use the Euler-Maruyama 
method for discretization of (0) with timestep At = 10~ 3 . To classify the asymptotic behaviour of the 
solution, we again define the “indicator” 


:=-Y 


9=1 


T/At 

E I 

fc=(T-l)/At 


r tf 


At 

Iv 


1/2 


(39) 
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Figure 3: Decadic logarithm of the indicator I a ,T, given by (J39J). for Monte Carlo simulations of the system 
([8]) with Q = 100 paths on the time interval [0,30] subject to the initial condition (1381) . with A = 1, (er, t) £ 
[0,2] x [0,2]. The dark blue regions (colour online) indicate “numerical flocking”. The solid line indicates 
the critical value t c given by formula o for A = 1 as a function of a. The number of individuals is: (a) 
N = 2; and (b) N = 20. 


where is the g-th path in the Monte Carlo simulation of (JT|) © at time tk = kAt. We say that numerical 
flocking takes place when I a T < 10” 2 . The background colour in Figure [3] encodes the decadic logarithm of 
the indicator, and the dark blue region indicates numerical flocking. We observe that the region of numerical 
flocking is only weakly influenced by the number of agents N. This is in agreement with the fact that 
the flocking condition C3 in Theorem [2] does not depend on N. The increased smoothness of the colour 
transition when IV = 20 is a consequence of the law of large numbers. For comparison, the critical value 
r c given by formula 1171) for A = 1 as a function of a is indicated by the solid line in both panels. The 
comparison with the numerical results suggests that the condition (1171) is far from optimal. 

4.4 Numerical study of the delayed Cucker-Smale system with multiplicative 
noise 

Finally, we present results of numerical simulations of the system CD 0 in the one-dimensional setting 
d = 1 with the communication rates i[>ij = if(\xi — Xj |) and %l> given by (CD- As in Section U~3l our goal is to 
characterize the asymptotic behaviour of the solutions in dependence on the parameter values, however, we 
are facing additional difficulties here. In particular, the asymptotic behaviour of the solution may depend 
nontrivially on the initial condition, as we show in Figure d Since a systematic study taking this effect 
into account is beyond the scope of this paper, we will impose the same type of initial condition for all our 
simulations. In particular, we prescribe constant zero value for the x-variables, 

Xi(t) = 0 for t £ (—r, 0], z = 1,2,..., iV. (40) 

For the v-variables we impose again the initial datum (1351) . 

We perform Monte Carlo simulations of the system CD 0, 0 with N £ {2,20} and (3 = 0.1 (strong 
coupling) and /3 = 1 (weak coupling). As in Section l4~3l we fix A = 1 and divide the domain [0,2] x [0,2] 
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Figure 4: Numerical simulations of the system (HJ) (J 2 ]) . j5]) with parameter values N = 2, X = 1, /3 = 0.1, 
cr = 0 and t = 1.75. Both simulations are performed on the time interval [0,120] with discrete timestep 
At = 10~ 3 . The initial condition for v is V\(t) = 1, V 2 (t) = —1 for t 6 (—r, 0] in both cases. The initial 
condition for x is: (a) X\(t) = —1, a ; 2 (t) = 1; (b) = 1, x 2 (t) = — 1. The plots show the velocities of the 

two agents (red and blue, colour online) as functions of time. 
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for (cr, r) into 50 x 50 equidistant (er, r)-pairs. For each pair of the parameter values we perform a Monte 
Carlo simulation with Q = 100 paths over the time interval [0, T] with T = 30. We use the Euler-Maruyama 
method for discretization of ® 0 with timestep At = 10 -3 . To classify the asymptotic behaviour of the 
solution, we again use the indicator (1TJ1) and say that numerical flocking takes place when I a>T < 10 -2 . 
The background colour in Figure [5] encodes the decadic logarithm of the indicator, and the dark blue region 
indicates numerical flocking. 

In the top left panel we indicate by an arrow the point (cr, r) = (0,1.75) that corresponds to the parameter 
setting in Figure 0 however, note that the initial conditions for Xi in Figure [4] differ from (l40l) . We see that 
the indicated point lies close to the boundary of the dark blue region, i.e., in the “transition zone” between 
numerical flocking and non-flocking. We hypothesize that this is why we were able to observe the two 
qualitatively different kinds of asymptotic behaviour in Figure [I] even if the initial datum for the v-variables 
is the same in both cases. Again, a systematic study of this hypothesis is beyond the scope of this paper. 

In Figure [5] we observe that the region of numerical flocking is only weakly influenced by the number of 
agents N. This is in agreement with the fact that the flocking condition (fTTl) in Theorem [2] does not depend 
on N. The increased smoothness of the colour transition when N = 20 is a consequence of the law of large 
numbers. On the other hand, we can distinguish two distinct types of patterns, one similar to Figure[2]for the 
strong coupling case j3 = 0.1 (Figures 0a) and 0b)), and a semicircular pattern for the weak coupling case 
/3 = 1 (Figures [51(c) and 0d)). In particular, the result for the weak coupling case is somewhat surprising 
- it suggests that for low levels of noise (cr < 0.6), introduction of intermediate delays (0.3 < r < 1.8) may 
facilitate flocking. This is further supported by Figure [6] where we plot sample solutions of ® (0, (0 for 
N = 2, /? = 1, cr = 0 (Figure 0a)), cr = 0.5 (Figure0b)) and three different values of the delay r £ {0,1, 2}. 
We observe that while for r = 0 and r = 2 the agents do not show tendency to converge to a common velocity 
during the indicated time interval, they exhibit numerical flocking for the intermediate value r = 1. We will 
call this observation time-delay induced flocking. 

Let us note that the results presented in Figure 0 do not contradict our analytical results. In particular, 
condition (fTTl) gives r < -\/2/4 = 0.35 if cr = 0 and r < (— 1/2 + yj ll/8)/4 = 0.17 if cr = 0.5, so it is only 
satisfied for the simulations in the panels corresponding to r = 0 in Figure0 Therefore, statement of Lemma 
[3] applies. The (expectation of) the Lyapunov function (1241) decreases in time for these two simulations. 

To gain a further understanding of the interesting phenomenon of time-delay induced flocking, we run 
systematic simulations of the system 0 0 .® with different values of f £ [0.5, 2.5], r £ [0,2] and cr £ 
{0,0.5}. We calculate the indicator Ig, T as in (f39l) with Q = 1 for a = 0 (there is no need to run more than 
one path for the case without noise) and Q = 100 Monte Carlo paths for cr = 0.5. The decadic logarithm of 
Ig, T is plotted in Figure [7] and we again use the threshold Ig tT < 10 -2 to define numerical flocking (dark blue 
regions in Figure [7]). We observe that there exists (for f sufficiently large) a region of intermediate values 
of r where numerical flocking takes place, while it does not for smaller or larger r values. Moreover, we see 
that noise has a disruptive influence on flocking (the dark blue region is smaller in Figure 0b) compared to 
Figure 0a)). 


5 Discussion 

We have studied a generalization of the Cucker-Smale model accounting for measurement errors through 
introduction of multiplicative white noise, and for delays in information processing. This has led to a system 
of stochastic delayed differential equations ® @. In Section 0 we have considered the communication rates 
between agents as given stochastic processes, and derived a sufficient condition for flocking , which we define 
as asymptotic convergence of the agents’ velocities towards a common value. The condition is given in terms 
of the critical delay that guarantees flocking as a function of the noise level. Our analysis is based on a 
construction of a suitable Lyapunov function for the system and a study of its decay. As a byproduct of the 
analysis, we obtain a sufficient condition for asymptotic convergence of delayed geometric Brownian motion. 
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Figure 5: Decadic logarithm of the indicator I a ,T, given by (EH), for Monte Carlo simulations of the system 
(jTj) (l2l) . © with Q = 100 paths on the time interval [0,30] subject to the initial conditions (1381) and (1401) . 
with A = 1 and (er, r) £ [0,2] x [0,2], The dark blue regions (colour online) indicate numerical flocking. The 
arrow in the top left panel indicates the point (er, r) = (0,1.75) that corresponds to the parameter setting in 
Figure 01 We use: (a) ft = 0.1, N = 2; (b) (S = 0.1, N = 20; (c) fj = 1, N = 2; and (d) /? = 1, N = 20. 
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Figure 6: Agents velocities v\(t),V 2 (t) in sample solutions of the system ([1} ([2]), © with N = 2, A = 1, 
/3 = 1 (weak coupling), on the time interval [0,30] subject to the initial conditions (l38l) and (l40l) . We use: 
(a) <t = 0; and (b) o = 0.5. 
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Figure 7: Decadic logarithm of the indicator Ig. T for simulations of the system (JTJ) ([2]), (O on the time 
interval [0,30] with A = 1, N = 2, (/3, r) S [0.5, 2.5] x [0,2]. We use: (a) ct = 0, Q = 1; and (b) a = 0.5, 
Q = 100. The dark blue regions (colour online) indicate numerical flocking. 


The second part of the paper is devoted to systematic numerical simulations. First, we perform Monte 
Carlo simulations of delayed geometric Brownian motion and evaluate its asymptotic behaviour based on a 
suitable “numerical indicator”. This led to the conclusion that the analytically derived sufficient condition 
for asymptotic convergence is qualitatively right - the convergence deteriorates with increasing noise level 
and delay. However, quantitatively it is far from optimal. Next, we simulate the Cucker-Smale type system 
with fixed communication rates and again compare with the analytical result. As before, the comparison 
shows that, while qualitatively correct, the analytical formula produces too restrictive critical delays. Finally, 
we simulate the full Cucker-Smale system with delays and multiplicative noise. We use two regimes for the 
dependence of the communication rates on the agents’ distances: the strong coupling regime, which leads to 
unconditional flocking in the “classical” Cucker-Smale model, and the weak coupling regime, where flocking 
may or may not take place. In the strong coupling regime the numerical picture is similar to the previous 
simulation with fixed communication rates. On the other hand, in the weak coupling regime we observe a 
somehow surprising behaviour of the system - namely, that an introduction of intermediate time delay may 
facilitate flocking. We call this phenomenon “delay induced flocking”. 

Our paper leaves several open questions. First of all, our analytical flocking condition is too restrictive 
compared to numerical results, so efforts should be made to improve it. Moreover, the analysis applied to 
the case when the communication rates are given and satisfying a certain structural assumption. This is 
in fact against the spirit of the original Cucker-Smale model where the communication rates depend on the 
mutual distances between agents. A possible extension of our analysis to this case remains an open problem. 
The main difficulty is due to the fact that it is not clear how to apply the classical bootstrapping argument 
that bounds the velocity fluctuations in terms of fluctuations in positions and vice versa. For the numerical 
part, it would be desirable to apply some multilevel Monte Carlo or importance sampling technique to obtain 
more accurate results. Moreover, the influence of the initial condition on the asymptotic behaviour should 
be studied. Finally, the interesting phenomenon of “delay induced flocking” deserves a detailed study, both 
from the analytical and numerical point of view. 
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